Absence of a Finite- Temperature Melting Transition in the Classical Two-Dimensional 

One-Component Plasma 



M. A. Moore 1 and A. Perez-Garrido 2 
Theoretical Physics Group, Department of Physics and Astronomy, The University of Manchester, MIS 9PL, UK 
2 Departamento de Fisica, Universidad de Murcia, Murcia 30.071, Spain 

Vortices in thin-film superconductors are often modelled as a system of particles interacting via a 
repulsive logarithmic potential. Arguments are presented to show that the hypothetical (Abrikosov) 
crystalline state for such particles is unstable at any finite temperature against proliferatio n of 
screened disclinations. The correlation length of crystalline order is predicted to grow as y/l/T 
as the temperature T is reduced to zero, in excellent agreement with our simulations of this two- 
dimensional system. 
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It has been commonly assumed for many years now 
that for the physically important case of particles mov- 
ing in two dimensions interacting with each other via 
a repulsive logarithmic potential (a situation sometimes 
called the two-dimensional one-component plasma prob- 
lem) one would have the usual phases expected on 
the KTHNY scenario. This scenario describes two- 
dimensional melting as a defect-mediated phenomenon 
(Halperin, Nelson M and Young 0) and is based on 
ideas of Kosterlitz and Thouless . It is supposed that 
the crystalline phase - a triangular lattice - melts at a 
continuous transition into an hexatic liquid due to the 
proliferation of dislocations. The hexatic liquid becomes 
an ordinary liquid at temperatures which permit the cre- 
ation of disclinations. This scenario is well-established 
for particles with short-range interactions Q] , but we will 
argue that particles interacting with a logarithmic po- 
tential behave completely differently. For them we can 
show that the crystalline state is unstable at any tem- 
perature against the proliferation of (screened) disclina- 
tions and as a consequence the system stays in the liquid 
state down to arbitrarily low temperatures. The ground 
state of the system is of course crystalline; the correla- 
tion length of short-range crystalline order is predicted 
to grow as ^/T/T as the temperature T approaches zero. 
Our numerical simulations reported here confirm this be- 
havior. 

The one-component plasma problem is of considerable 
physical significance as it relates to the thermodynamics 
of vortices -"the particles" - in thin film superconduc- 
tors. For thin enough films the screening length in the 
intervortex potential may be greater than the transverse 
dimensions of the film, which makes the logarithmic po- 
tential an accurate approximation for the potential. Most 
papers on thin film superconductors assume the vortices 
have a freezing transition at low enough temperatures 
(for a review see Ref. ||), although clear experimental 
evidence for this is lacking. For a contrary view, how- 
ever, see H and references cited therein. 

For particles interacting with a repulsive potential a 
device is needed to stop them escaping to infinity. In 
numerical studies of two-dimensional melting the most 
commonly used device is periodic boundary conditions. 



Unfortunately the use of this boundary condition with 
either short-range interactions or with the logarithmic 
interaction Q produces an apparently first order tran- 
sition between the crystal and liquid states rather than 
the KTHNY scenario. (For a review of early work on 
short-range interactions see for some more recent 
work see Refs. Q or (it))). This is probably a finite size 
effect: studies on systems with over 60,000 particles in- 
dicate that the van der Waals loops associated with the 
apparent first-order transition shrink in these very large 
systems ||, [[To). We ourselves have found that placing 
the particles on the surface of a sphere is very effective 
for short-range interactions |ll]] : no van der Waals loops 
occur with this topology and the results obtained even 
with modest numbers of particles are in excellent agree- 
ment with expectations based on KTHNY theory. As a 
consequence the numerical work which we are reporting 
in this paper has been carried out for the two-dimensional 
system represented by the surface of a sphere. 

The ground state configuration of the particles on the 
sphere has to contain at least 12 disclinations (5-fold 
rings) by Euler's theorem. We have made extensive stud- 
ies of these ground states and discovered that for larger 
systems the disclinations are screened by lines of disloca- 
tions (|], K3]. These defects within the crystalline state 
seem to overcome the problem of the spurious first or- 
der transition induced by finite size effects when periodic 
boundary conditions are employed and so enable one to 
get results closer to those obtaining in the thermody- 
namic limit. It is noteworthy that an early simulation of 
the one-component plasma on the surface of a sphere |l~3| ] 
did not find a finite temperature phase transition either, 
in agreement with our results. 

The Hamiltonian for particles moving on the surface 
of the sphere interacting via a logarithmic potential is 
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where is the position of the ith particle on the sur- 
face of the sphere, R is the radius of the sphere and J 
is a measure of the magnitude of the repulsive forces be- 
tween the particles. The key feature which distinguishes 
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the logarithmic potential from other potentials is that all 
the stationary states of H have zero dipole moment jl4| . 
This was proved by noting that the force on the ith par- 
ticle due to all the others must be directed radially for 
any equilibrium configuration since otherwise the particle 
would move along the sphere, so 
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By multiplying both sides by r.; one can show that /$ = 
(N — 1)/2R 2 where N is the total number of particles. 
By summing Eq. (Q) over all i and using the fact that 
Vi — Tj is antisymmetric in i and j, it follows that the 
dipole moment, X)i r ij i s zer O- No other potential has 
this feature and it has important consequences. 




FIG. 1. Dislocation series screening a five-fold disclination. 

Our basic contention is that thermally excited screened 
disclinations will destroy crystalline order for particles 
interacting with each other via a logarithmic potential. 
The energy cost of an unscreened disclination is O(N) 
H and such a disclination will not be thermally cre- 
ated in a crystalline state. However, disclinations can 
be "screened" by a cloud of dislocations and it turns out 
that the energy cost of such screened disclinations can 
be much smaller, of 0(ln(N)) for non-logarithmic po- 
tentials and 0(1/N) for the logarithmic potential. The 
phenomenon of screening of a disclination by dislocations 
is well-known Q, |J, and is illustrated in Fig. 1. The 
figure shows a five-fold coordinated ring - a disclination 
- screened by five lines of dislocations where the dislo- 
cations are spaced a distance cIq apart; lo is the lattice 
spacing. The strain field of the central disclination can 
be largely cancelled by that arising from the lines of dis- 
locations. The resulting strain field from the dislocations 
along a line may be approximated at large distances by 



that which arises from a positive and negative disclina- 
tion at each end of the line with a "disclination" charge 
of size q s = 1/c (where q s — +1 for a fivefold disclina- 
tion) . If we consider a line of dislocations with c = 5 and 
five lines as in Fig. 1 then the strain field of the central 
disclination is exactly screened away. As shown in Ref. 
H the contribution of the disclinations at the other end 
of the lines can be made arbitrarily small by allowing 
the spacing of the dislocations to increase with distance 
r from the disclination as clr) = 5 + S/log(r/S) with 
the condition that g(0) — 0; S is the size or scale of the 
screened disclination. Then the residual charge associ- 
ated with the screened disclination can be made as small 
as 0{Iq/S) but, as we shall show, must be as small as 
O((lo/ S) 2 ) for the special case of the logarithmic poten- 
tial. 

First let us review some features of two-dimensional 
continuum elasticity theory p|. Small strains ity(r) 
are related to the stress field by Hooke's law, ov,- = 
B5ijUkk+2ii(uij—5ijUkk/2), where B is the bulk modulus 
and [i is the shear modulus. In the presence of topological 
defects it is convenient to introduce the Airy stress func- 
tion, X: defined by <jjj = en-CjidkdiX- A fivefold discli- 
nation is defined by a change in bond angle 2it/6 when 
a path encircles the defect. Dislocations are defined by 
their Burgers vector density field b(r) which for the dis- 
locations in Fig. 1 points perpendicular to the line upon 
which they lie. The stress field is related to the densities 
of disclination s(r) and dislocations via 
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eifeV fe 6i(r) = s(r), 
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where Yi = 4B/i/(B + fi). For a single disclination at 
the origin as in Fig. 1, s(r) = (27r/6)<5(r). s(r) can 
be regarded as a total disclination density made up of 
a "free" disclination density s(r) and a "polarization" 
contribution — ejfcVfc&j from dislocations. The energy of 
the screened disclination expressed in Fourier space is 
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The expectation || for non-logarithmic interaction po- 
tentials is that the screening can at best make the Fourier 
transform of s, s(q) = qlof(qS) when the amount of 
"disclination charge" within a region of radius S around 
the disclination is of 0(Iq/S). Substituting this form for 
s(q) into Eq. one finds that in a system of N par- 
ticles the energy of the screened disclination would be 
of 0{lnN) - which explains why screened disclinations 
would be unlikely to modify the KTHNY scenario for 
non-logarithmic interactions. 

The change in the particle density due to the presence 
of the topological defects is, when Fourier transformed, 
given by 



Ap(q) = -7m»(q) = ^<2 2 x(q), 



(5) 
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where n is the number density of the particles. A fea- 
ture of the logarithmic potential is that for it the bulk 
modulus B is not constant pa but diverges at small 
wavevector: B(q) — 2itJn 2 /q\ The shear modulus is 
well-behaved: p, = Jn/8. Equations (Q) and (||) can still 
be used if the replacements Y% — > 4/i and B — > B(q) are 
employed. At small wave- vector it follows that for the 
logarithmic potential 
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Only for the logarithmic case does a finite small q 
limit exist for the density change Ap associated with a 
screened disclination. 

We now exploit the fact that all stationary states of 
H have vanishing dipole moment to show that for the 
case of logarithmic interactions between the particles the 
screening of the disclination is more efficient than for 
non-logarithmic interactions. The Fourier transform of 
the particle density is defined by 
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(Formally, as our system is the surface of a sphere rather 
than a plane we should use spherical harmonics rather 
than plane waves, as was done in Ref. [n, but the dis- 
tinction is unimportant for our argument). Consider now 
a state which differs from the groundstate by the pres- 
ence of a screened disclination of size S. The density 
difference Ap(q) of the two states must differ as q — > 
by terms of 0(q 2 ) ; (if one of the states had had a dipole 
moment then one can see from expanding the exponential 
in Eq. (0) that Ap would have been of 0(q)). Eq. @ 
implies that s(q) = q 2 lof{qS). By Fourier transforming 
s(q) one can then show that the "disclination charge" 
within a distance S of the center of the disclination is of 

Substituting this form for s into Eq. (|4j) one finds 
that the energy of the screened disclination is of order 
J{Iq/S) 2 . By increasing the scale S it can be made ar- 
bitrarily small. At a temperature T a region of linear 
extent £, where J{lo/£,) 2 — T, will be unlikely to contain 
a disclination and so £ will be a measure of the short- 
range crystalline order present in the system at temper- 
ature T. £ diverges as \JTJT as T — > 0. This means that 
by investigating numerically the structure factor one can 
find from the widths of its peaks the correlation length £ 
and its temperature dependence will tell us whether the 
arguments above are valid. 

We studied using molecular dynamics, specifically a 
velocity Verlet algorithm |lq] , a system of N particles 
confined to move on the surface of a sphere and interact- 
ing with the logarithmic potential. Reduced units were 
used, i.e. m — k^ = R = J = \, where m is the mass 
of the particle and k& is the Boltzmann constant. The 
acceleration an of the ith particle equals fi/m, where fi 
is the force produced by the other particles on the zth 



particle. After a small time interval St the position of 
the particle will be x, = rj(i) + V{(t)St+ ^&i{t)5t 2 , where 
Vi(t) is the velocity of the particle. In general Xj will 
not lie on the surface of the sphere. The zth particle is 
brought back to the surface by acting on it with a ficti- 
tious force — 2XiTi(t) where 
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ri(t)-xi- J[rj (i)-xi] -R 2 
BJSt 2 



R 2 
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Then, the velocity Verlet algorithm updates particle po- 
sitions using the equation: 



Ti(t + St) =x i -X i r i (t)6t 2 . 



(9) 



We chose St = 0.005(mi? 2 /^) 1/2 - The velocities of the 
particles were chosen from a Boltzmann distribution ap- 
propriate to the temperature T and were re-selected at 
equally spaced time intervals Jl7| . The system was equi- 
librated at high temperatures, then the temperature was 
slowly reduced. For each temperature we determined the 
structure factor which is related to the Fourier transform 
of the pair correlation function h(r) by ]lq|: 



S(q) = 1 + 2-KpRf \ h(R9) sm6J (qR9)d9, (10) 
Jo 

where Jq is the Bessel function of zeroth order. This 
adaption of the conventional relation to particles mov- 
ing on the surface of the sphere is valid provided q is 
of O(l) and not of 0{\/R). Peaks in the structure fac- 
tor grow as the temperature is reduced at wavevectors q 
corresponding to the reciprocal lattice vectors |G| of the 
triangular lattice expected for the groundstate. The cor- 
relation length, £, is the inverse of the width of the first 
peak of the structure factor. To determine it, we fitted 
the first peak to a Lorentzian curve. 
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FIG. 2. Correlation length as a function of T~ 1/2 for 1442 
particles (solid circles) and for 2252 particles (empty dia- 
monds) for the logarithmic potential. 
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We studied systems of 1442 and 2252 particles. The 
simulation time for each temperature was 100,00(Wt. In 
Fig. [2] £ is plotted against y/l/T. The vertical line is 
drawn where other authors found a first-order melting 
transition using periodic boundary conditions JtJ. The 
predicted behavior £ oc y/T/T is clearly seen in Fig. |^. 
When the temperature was reduced to T — 0.01 the cor- 
relation length for the system of 1442 particles reached 
the system size and stopped growing as the temperature 
was reduced further. However, the simulation with 2252 
particles indicates that this levelling off is just a finite 
size effect. True equilibration in numerical studies of two- 
dimensional melting phenomena is always problematic ||] 
and may be the cause of the scatter in Fig. |[ 



p c = 1, a value obtained by other authors [|19|;|9| working 
at the temperature which we used, T = 1. The slope 
of the straight line is —1 which corresponds well with 
KTHNY predictions. Note again that finite size effects 
cut off the growth of £ when it is of O(R). Thus sim- 
ulations on the sphere do produce for a non-logarithmic 
potential the expected crystalline phase. 

In summary, we have shown that thermal excitation of 
screened disclinations removes at non-zero temperature 
the crystalline phase of the vortex system. Numerical 
simulations have confirmed our prediction that as the 
temperature is lowered the correlation length of short- 
range crystalline order should grow as £ oc y/l/T. 
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FIG. 3. Log-log plot of ln£ versus (1/p 6 - l) - 36963 for 5882 
particles interacting with a 1/r 12 repulsive potential. The 
slope of the straight line is —1 according to KTHNY expec- 
tations. 

The apparent absence of a finite temperature phase 
transition for particles interacting with a logarithmic po- 
tential cannot be attributed to the fact that we have done 
the simulation for the two-dimensional geometry repre- 
sented by the surface of a sphere. We can demonstrate 
this by simulating particles interacting with a 1/r 12 po- 
tential also moving on the surface of a sphere. A finite 
temperature melting transition of KTHNY character is 
seen. We had already found indications of such a melt- 
ing transition flT|| but the numbers of particles used in 
that reference were rather small. Using the Verlet al- 
gorithm described above we are able to simulate much 
larger systems eg. 5882 particles. The calculations for 
the short-range potentials run faster than with the log- 
arithmic potential as it is possible to use look-up tables 
of nearest neighbors. In the KTHNY picture, the corre- 
lation length has the following density dependence along 
an isotherm for a 1/r 12 potential p9[: 
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